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Abstract. Two years of neutral measurements by IBEX-Lo have yielded several direct observations 
of interstellar neutral helium and oxygen during preferred viewing seasons. Besides the interstellar 
signal, there are indications of the presence of secondary neutral helium and oxygen created in 
the heliosphere. Detailed modeling of these particle species is necessary to connect the measured 
fluxes to the pristine local interstellar medium while accounting for loss and production of neutral 
particles during their path through the heliosphere. In this contribution, global heliosphere models 
are coupled to analytic calculations of neutral trajectories to obtain detailed estimates of the neutral 
distribution function of primary interstellar helium atoms in the heliosphere, in particular in the 
inner heliosphere. 
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1. INTRODUCTION 

The neutral flux maps accumulated by NASA's Interstellar Boundary Explorer (IBEX) 
since the start of its science operations in early 2009 contain detections of direct in- 
terstellar neutral helium and oxygen atoms [1]. These detections are centered around 
time periods (in Spring) when IBEX scans in the direction of the difference vector be- 
tween Earth's orbital velocity and the gravitationally deflected and accelerated particle 
velocity. The main particle streams are consistent with a far-away velocity (in the local 
interstellar medium, LISM) determined by [2] who obtained the helium flow field by 
careful analysis of the Ulysses/GAS measurements (see also [3] for a synopsis of the 
values from different techniques). Further analysis of the IBEX data reveals the pres- 
ence of secondary oxygen and secondary helium in the measurements. Those particles 
are created through charge exchange of He + , He ++ , or + with heliospheric neutrals 
(H, He, O) and hence are in a different velocity state, with different neutral trajectories. 

Calculating primary and secondary neutral fluxes at IBEX through modeling requires 
combining a global heliosphere model treating the dominant neutral H and H + (pro- 
tons) self-consistently to obtain a realistic heliospheric plasma background. On this 
background, heavy neutrals like He and O can be simulated with Monte-Carlo- style 
kinetic codes [e.g., 4] where primary interstellar neutrals are injected at an outer bound- 
ary, sampled from a thermal Maxwellian, and transported through the heliosphere while 
subjected to loss processes (interaction with the background plasma and solar photons). 
Secondary neutrals are generated by charge exchange everywhere in the heliosphere, 
and after creation they are transported as well, while subject to losses. 

This paper pursues an alternative method for calculating heavy neutrals, using Kep- 
lerian equations to describe the entire trajectory while conserving quantities like energy 



and angular momentum to maximize the accuracy of the results . Using conserved quanti- 
ties connects a 1 AU location (IBEX) with the pristine ISM (LISM) in one algebraic step. 
With this efficient method, the goal then is to carry out precision calculations, where the 
entire phase space distribution function at a desired point is calculated as accurately as 
desired without the steep computational cost that a Monte Carlo method typically incurs 
for a similar accuracy. 

2. CONSERVED QUANTITIES ALONG TRAJECTORY 

The Sun exerts a central force on neutral particles in the heliosphere, and the neutral 
trajectories hence are Keplerian orbits with the Sun at the focus. The motion proceeds 
in a central potential -GMq/t. More generally, the potential constant is defined as 
ffx = GMq(\ - ju), with G the gravity constant, Mq the solar mass (GMq = 887.6 
AU km 2 s~ 2 is used here), and ji the ratio of radiation pressure force to gravitational 
force. While /i 1 for hydrogen, /i = for heavy neutrals, with solar gravity acting 
unopposed. 

The motion of an individual neutral particle is confined to a plane (orbital plane), 
containing the Sun, the particle, and its velocity vector. There are several conserved 
quantities along the entire trajectory. The most important ones are the following seven, 
which are not all independent from each other: Angular momentum, Laplace-Runge- 
Lenz vector, and total energy. As those all involve the particle mass or its square as 
prefactors, the latter can be made part of the conservation constants, to arrive at the 
following set of conserved quantities: 

— * 

I = rxv (specific angular momentum) (1) 

— * 

a = vxl — f^r (eccentricity vector) (2) 

V 2 fa 

E t ot = (total specific energy) (3) 

2 r 

With an appropriate coordinate system, two angular momentum components and one 
eccentricity vector component can be made to identically vanish everywhere. The ec- 
centricity vector a points in the direction of perihelion, and its magnitude is related to 
the orbital eccentricity e by \a\ = ef^. 

The conserved quantities can be used in at least two ways: In a direct way, when a par- 
ticle's location and velocity is specified at one time. Then, all conserved quantities are 
immediately determined, and the entire future (forward) trajectory and past (reverse) tra- 
jectory can be calculated with formulae from celestial mechanics. With the knowledge of 
the conserved quantities, any desired radial distance, azimuthal angle ("true anomaly"), 
point in time, or radial velocity, for example, can be used as input to determine all other 
variables accurately, with the conserved quantities conserved by design. A natural ap- 
plication of this is backtracking a particle to a large distance where it is sure to be in the 
LISM. In this way, a particle measured by IBEX with a certain direction and energy can 
be backtracked to where it came from, and the travel time from a reference distance to 
IBEX can be calculated algebraically as well. 



In the second way, the mixed problem of knowing a particle's position at one point in 
time, and a particle's velocity vector at infinite distance where it started out, is solvable 
as well. Here, two sets of conservation equations (l)-(3) are established, and treated as 
a system of equations that, when solved, yields the missing variables. An application of 
this problem has been published quite some time ago in the context of characterizing 
the velocity state of neutrals in the inner heliosphere when a cold interstellar neutral 
distribution (known velocity at infinity) is assumed [5]. Equation 48 of [5] gives the 
algebraic formula for the angular momentum and hence velocity vector at the location 
of interest. Each location is reached by two neutral trajectories (two solutions), a direct 
path and a longer indirect path. Typically, the indirect path takes the particle closer to 
the Sun, and for a longer time, and is hence more affected by losses than the direct 
path; see Figure 2a below for an illustration which also contains the respective values 
for eccentricity and specific angular momentum (the latter in units of AU km/s). 

3. GRAVITATIONALLY DEFLECTED PHASE SPACE DENSITY 

In the pristine LISM a Maxwellian phase space distribution (PSD) can be assumed for 
the interstellar neutral helium atoms. Through the action of gravity, the PSD shape at a 
point in the inner heliosphere is different from Maxwellian, even if charge-exchange col- 
lisions are ignored. The PSD maximum is determined by using the bulk speed vector of 
the LISM [2] as input and solving the mixed problem as outlined above. The entire PSD 
is obtained by scanning the velocity space around this maximum in all velocity space 
directions; for each new velocity, this "direct problem" (see above) is solved, backtrack- 
ing the particle to yield the velocity at infinity. By virtue of Liouville's theorem, the 
PSD of the Maxwellian for this velocity equals the one of the phase space point under 
investigation. For practical purposes, such scanning is ended when the PSD value falls 
below 1(T 4 of the maximum value, which is on the order of 4 thermal velocities away 
from the bulk speed in the LISM. 

Consider an example of this procedure, with a sample point of interest chosen 2 AU 
away from the Sun, in a downwind- sidewind direction that is not yet in the helium 
focusing cone. The location coordinates r=(- 1 .0, 1 .73, 0.0) in a frame of reference where 
the x axis points into the LISM flow, combined with the LISM bulk velocity (-26.3, 0.0, 
0.0) km/s [2], give the direct-path velocity solution (-36.0, -16.8, 0.0) km/s and the 
indirect solution (-3.4, +39.6, 0.0) km/s. The situation is illustrated in Figure 2a. In the 
entire velocity space at r, there are two maxima of value 1 at these two velocities, each 
surrounded by a finite volume where the normalized phase space density remains larger 
than 0.001, an arbitrary value chosen to envelop the most important pieces of phase 
space. Figure 1 shows perpendicular cuts through these two 3D velocity space objects, 
all taken at planes that contain the respective point of maximum. The shape of these 
cuts suggests that the deflected direct PSD has a form similar to a thick contact lens; 
the indirect PSD is more irregular yet confined to a smaller volume. For comparison, 
the LISM Maxwellian contours are concentric circles, with the 0.001 contour having a 
radius of 14 km/s. 

The procedure so far disregards any losses interstellar helium might experience on 
its way from the pristine LISM to a point in the inner heliosphere. Photoionization 




FIGURE 1. Phase space densities of primary neutral interstellar helium at the location (-1.0, 1.7) 
slightly downwind of the Sun and away from the symmetry axis (2 AU distance). The PSD is a 3D 
object; the left column presents v x -v y cuts through the center (v z =0 in both panels; note the v x axis shift 
between the panels). Analogously, v y -v z cuts taken at the centers v x =-3 km/s (indirect PSD) and at ^=-36 
km/s for the direct PSD constitute the right column of the figure. The maximum value is normalized to 1 
(no losses included); the contour levels are spaced logarithmically with 3 lines per dex. 



is the dominant loss; its local rate is proportional to the inverse square of distance to 
the Sun. Charge exchange (ex.) channels are the next most important losses; their rate 
is proportional to the background plasma density, the relative plasma-neutral speed, 
and ex. cross section. The ex. channels are loss of helium due to encounters with 
plasma protons, with He + ions in the outer heliosheath, and with He ++ ions (alpha 
particles) in the solar wind. Some reasonable abundances [6] are chosen to represent 
helium ion densities and velocities with the help of the proton plasma background of 
an axisymmetric global heliosphere model. The introduction of the numerical plasma 
background presents an increase in the amount of calculations: Now, the orbital position 
and velocity at each background grid point need to be calculated to evaluate the local 
loss rate. With the conserved quantities available, this still is a one-step calculation; it 
is just the number of such calculations that goes up, roughly to the order of one linear 
dimension of the background plasma grid. 

In keeping with the above example, Figure 2 shows the modified v x -v y cuts of the PSD 
when all four loss processes are taken into account. The direct PSD is lower now, with 
a new maximum of 0.61 (the first contoured level is 0.46). The PSD has a smaller slope 
now, so that the 0.001 contour location is not very different from before. In contrast to 




this, the indirect path takes the particle closer to the Sun, and the corresponding PSD 
region suffers more dramatic losses, with even the 0.001 contour surrounding a much 
smaller volume. The maximum achieved at the center is 0.14, indicative of strong loss. 

With the PSD characterized in such a detailed fashion, it is easy to integrate moments 
to arrive at primary neutral density, bulk velocity, and temperature, also separate for 
direct and indirect populations if desired. The PSD in the wake of the Sun, where a 
focusing cone exists, has a complex shape, degenerating into a ring on the downwind 
symmetry axis. There, higher moments do not characterize the actual neutral PSD well. 

To complete the picture and make detailed connections with observations, it will be 
necessary to add secondary neutral particles to the calculations. They are the neutral 
product of charge exchange of a helium ion with a neutral H atom (supplied by a 
background global heliosphere model) or a neutral He atom (as calculated above). He 
ions are part of the solar outflow and the ISM inflow, as mentioned above for the ex. 
loss processes of primaries. Naturally, the pickup ions generated in primary ex. events 
are a further source of He ions available for the production of secondary neutral He. 
It is often assumed that the level of pickup ions is negligible when contrasted with the 
solar wind and interstellar helium ion densities, but this should be verified in detail for 
typical global heliospheric models. In case the pickup ions contribute significantly to 
the secondary helium, techniques for their inclusion in the global models will include 
location-dependent extra source terms for the helium ion fluid. 

While the source of primary helium is at infinity, the source of secondary helium lead- 
ing to a particular phase space point is distributed along the entire trajectory. Therefore, 
the secondary PSD will involve by necessity an integral along the trajectory, with pri- 
mary He calculations at each point of the trajectory. In this sense, non-primary particle 
calculations are a hierarchical process where "higher orders" are derived from "lower 
orders", and therefore call for more involved calculations. These issues will be treated 
in future papers. 



4. CONCLUSIONS 



Using the Kepler equations of celestial mechanics is advantageous when calculating 
trajectories in the heliosphere of neutral heavy atoms originating in the local interstellar 
medium. The associated conserved quantities enable efficient, one-step calculations of 
particle locations and velocities that are accurate by design. 

One example was presented to show the accurate output of the method, emphasize 
the role of gravitational deflection on the helium phase space distribution, and discuss 
the changes that further result from loss processes. There are two maxima in the helium 
PSD in the heliosphere, but typically the indirect path contribution suffers more severe 
losses. 

The primary neutral helium PSD can be calculated in this way at any arbitrary point in 
the heliosphere, and its various moments will yield maps of effective temperature, bulk 
velocity, and density. The latter will be elevated in the region of the helium focusing 
cone, which is also where direct and indirect solutions start to merge together and 
therefore become equally important. The general helium filtration by the heliospheric 
interface, by solar radiation and the inner solar wind is more pronounced than earlier 
studies suggested (e.g., [4]) in which only charge exchange with protons was taken into 
account as loss process. These issues will be studied in more detail elsewhere. 
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